R5R Tutorial

Introduction

In this tutorial you will be introduced to r5r: an interface to the R5 routing engine developed by Conveyal. R5 allows you to do “rapid, realistic” routing on multimodal networks. It’s very flexible, enabling intermodal analyses and has the key strength of retaining the full detail of the public transport timetable. If you’re more familiar with python, you can consider using r5py instead - although it may not have all of the same functionality.

The tutorial consists of the following main steps:

1. Installing r5r
2. GTFS preparation
3. Building a r5r model
4. Preparing origin-destination data
5. Calculating and exporting a TTM

In addition to r5r, you will be exposed to several other useful R packages: dplyr (for working with data tables), sf (for working with geodata), and tidytransit (for working with GTFS data). These are all very powerful packages that can make your life a lot easier if you learn how to use them! You’re encouraged to look into the documentation of these packages for more information.

Tip: you can modify this file and use it as a template to help you complete your project ;)

Installing r5r

All the packages listed below (including r5r) can be installed by clicking “packages > install” on the sidebar and searching for the package name or by executing the following command: install.packages("package_name_here")

Since R5 (the routing engine itself) is written in java, we need to install the Java Development Kit (JDK) to get it to work. Specifically we need version 21. You can install it here: https://www.oracle.com/de/java/technologies/downloads/#java21

Once that’s done, we can load all the packages we’ll be using today.

library(tidyverse)
library(sf)
library(tidytransit)

options(java.parameters = '-Xmx4G') # we need to allocate RAM for r5r to work. You can increase this depending on your system resources.
library(r5r)

For the code to run correctly, you need to set the “working directory”. This should be the repository (folder) you downloaded / cloned from github.

setwd("C:/you_need_to_update_this_path/ilutm-r5r-tutorial/")

Error in setwd("C:/you_need_to_update_this_path/ilutm-r5r-tutorial/") : 
  cannot change working directory

GTFS preparation

• Download data from: https://gtfs.de/ - pick “Deutschland gesamt”
• Move the zip file into the /raw_data/ folder
• Using the full feed for all of Germany will require a lot of resources. Since we’ll be analyzing a single city - lets filter the feed so that it is smaller. We’ll be using tidytransit for this.
• Let’s first load in geodata representing the area we’ll be trimming to (in this case - Munich). We’ll use st_read() from the sf package for this. It can handle all typical geospatial file types (.shp,.gpkg,etc.)

area <- st_read("./raw_data/munich_admin.gpkg")

Reading layer `boundary_administrative_munich' from data source 
  `C:\Users\barte\Documents\R\ilutm-r5r-tutorial\raw_data\munich_admin.gpkg' using driver `GPKG'
Simple feature collection with 25 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 11.36078 ymin: 48.06162 xmax: 11.72291 ymax: 48.24812
Geodetic CRS:  WGS 84

ggplot()+
  geom_sf(data  = area)

• It’s good practice to buffer our study area to avoid any edge effects. Let’s buffer our area by 5 km. NOTE: in this tutorial we’re just going to be looking at accessibility within Munich. If your destinations extend beyond the border of your study area, you may need to buffer your study area significantly! (>>5 km)
• We’ve seen how to do this in QGIS, let’s use sf this time

area_buffered <- st_transform(area,crs = 25832) # switch crs to projected (ETRS89 / UTM zone 32N)

area_buffered <- st_union(area_buffered) # join geometries 

area_buffered <- st_buffer(area_buffered, 5000) # buffer by 5km

area_buffered <- st_transform(area_buffered,crs = 4326) #switch back to WGS84 for filtering the GTFS feed


ggplot()+
  geom_sf(data = area_buffered)+
  geom_sf(data = area)

• Let’s load in the GTFS feed we downloaded and trim it.

gtfs_raw <- read_gtfs("./raw_data/latest.zip",encoding = "UTF-8") # path should point to the GTFS feed you downloaded

|--------------------------------------------------|
|==================================================|

gtfs_filtered <- filter_feed_by_area(gtfs_raw,st_bbox(area_buffered))

• We can see that the number of stops in the feed has been significantly reduced.

• Let’s save the filtered feed. R5 will use this to build a network, you can also use it for other public transport analyses (refer to previous tutorials). Tip: tidytransit has a lot of great functionality for this!

write_gtfs(gtfs_filtered,"./r5r_model/gtfs.zip")

Building a r5r model

• Before running any more code we need to supply R5 the data it needs to build the network.
• So far we have prepared the GTFS data, but we also need a street network.
• R5 uses an extract of OSM data to build the street network. To do this, it needs to be provided with a .pbf file. This can be found on Geofabrik or https://slice.openstreetmap.us/.
• OSM by the slice is recommended because it has a cool name and because you can select the exact area you’re interested in (reducing network size, improving performance)
• Just as when we were trimming the GTFS feed, extract an area larger than the study area

• Save the .pbf file to the /r5r_model/ folder. The folder should look like this:

• We can now tell r5 to build the network.

• Tip: if you are sure you aren’t modifying the GTFS or .pbf data any further, you can turn overwrite = F to avoid rebuilding the model (faster)

r5r_core <- setup_r5(data_path = "./r5r_model", verbose = FALSE,overwrite = T)

Preparing origin-destination data

• In this tutorial we’ll use public transport stops as our origins and destinations.
• As we’ve seen in a previous tutorial, we can get this information directly from the GTFS feed. Let’s do it in R instead of QGIS this time
• As we discussed in lecture, a single “stop” may be represented by multiple in a GTFS feed (e.g., a bus stop on both sides of the street).
• We can take some steps to simplify this. We will group all stops that share the same value for parent_station and then keep a single coordinate by taking the average. Be aware that solving this in a line of code is often too good to be true. GTFS data often contains mistakes, and may require more careful cleaning!

stops <- summarize(group_by(gtfs_filtered$stops, parent_station),
                   stop_lat = mean(stop_lat),
                   stop_lon = mean(stop_lon),
                   parent_station_name = first(stop_name))  
stops

• let’s convert the table to geodata and filter it down to our study area.

stops <- st_as_sf(stops,coords = c("stop_lon","stop_lat"),crs = 4326,remove = F)

stops <-st_filter(stops,area) #apply spatial filter

stops

Simple feature collection with 1125 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.38896 ymin: 48.06913 xmax: 11.71527 ymax: 48.22
Geodetic CRS:  WGS 84

• Finally, let’s format our data. The origins and destinations should have an id column. Let’s add a new column called id, and set the values to equal parent_station.

stops <- st_as_sf(stops,coords = c("stop_lon","stop_lat"),crs = 4326,remove = F)

stops <-st_filter(stops,area) #apply spatial filter

stops

Simple feature collection with 1125 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.38896 ymin: 48.06913 xmax: 11.71527 ymax: 48.22
Geodetic CRS:  WGS 84

Calculating and exporting a TTM

The main functionality of r5r is the ability to calculate a travel time matrix. Knowing the travel time between pairs of points is useful on its own. It can also be used if you want more control over your accessibility calculations. For example, if you calculate a travel time matrix between pairs of public transport stops, you can calculate the access (origin –> first public transport stop) and egress (last public transport stop –> destination) legs of the trip independently. In other words, combining the output of r5r with the QGIS-based accessibility calculations we did in earlier tutorials.

Let’s give it a shot:

tip: there are many other settings you can change, you’re encouraged to read the documentation of r5r to see what is possible! The documentation has many nice examples!

• stops: our origins are the public transport stops
• destinations: our destinations are the public transport stops
• mode: the primary mode is public transport (which includes walking for access and transfers)
• mode_egress: the egress mode is walking
• departure_datetime: start time
• time_window: A travel time matrix will be calculated for each minute, for time_window minutes after the departure time. This allows for us to account for the variability in travel times. Imagine we ran our analysis only at 8:00 and we just missed a departure of an S-Bahn at 7:59, the travel time that we’d calculate would be severely penalized by an excessive wait time. Calculating travel times over a period of time lets us mitigate this. This parameter can be reduced (default is 10) as it can be computationally expensive. It’s increasingly important in suburban and rural areas where service frequencies are low.
• percentiles: this corresponds to the time_window. a value of 50 will report the median travel time. We can specify multiple percentiles by using a vector c()
• max_walk_time: max time for access, egress, and transfers by walking
• walk_speed: (kph)
• max_ride: maximum # of transfers

TTM<-
travel_time_matrix(
  r5r_core,
  origins = stops, 
  destinations = stops, 
  mode = c("TRANSIT"),  
  mode_egress = "WALK", 
  departure_datetime = as.POSIXct("26-05-2025 9:00:00", format = "%d-%m-%Y %H:%M:%S"), #its important that the date and time is within your GTFS feed
  time_window = 60L, # 
  percentiles = c(10,50),
  max_walk_time = 10,
  max_trip_duration = 30L,
  walk_speed = 4,
  max_rides = 3,
)

TTM

The impact of the time_window parameter can be quite significant on an analysis! In this example we calculated:

• the 10th percentile travel time (meaning 10% of trips are at least this fast)
• the 50th percentile travel time (meaning 50% of trips are at least this fast)

The proportion of OD pairs that had a 50th percentile travel time > 30 min and a 10th percentile travel time <= 30 min:

TTM <- mutate(TTM,tt_diff = travel_time_p50 - travel_time_p10)

print(round(sum(is.na(TTM$tt_diff)) / length(TTM$tt_diff),1))

[1] 0.4

For those with a 10th and 50th percentile travel time <= 30 min, the distribution of the difference between the two is as follows:

TTM$tt_diff%>%na.omit()%>%summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.000   4.000   3.675   4.000  17.000 

TTM<-TTM%>%select(-tt_diff)

Let’s get the travel time matrix ready for an export. We will restructure the data so that the geometry of the destination stops (to_id) is added to each row. We will also add the name of the origin station so that the data is easier to work with. This is a similar operation to a VLOOKUP in excel. We use the from_id and to_id fields as they represent the parent_station.

For more guidance on dplyr functions (like left_join()) check out this cheat sheet

# join the geometry of the destination station
geometry_to_join <- select(stops,id) # drop unnecessary cols
geometry_to_join <- rename(geometry_to_join,to_id = "id") # rename "id" to "to_id" so that the data is joined as expected

TTM <-
  left_join(TTM, 
            geometry_to_join,
            by = "to_id")


# join the station name of the origin station

names_to_join <- select(stops,parent_station_name,id) # drop unnecessary cols
names_to_join <- st_drop_geometry(names_to_join) # drop geometry
names_to_join <- rename(names_to_join,from_id = "id") # rename "id" to "from_id" so that the data is joined as expected


TTM <-
  left_join(TTM, 
            names_to_join,
            by = "from_id")

TTM <- st_as_sf(TTM)

TTM

Simple feature collection with 295492 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.38896 ymin: 48.06913 xmax: 11.71527 ymax: 48.22
Geodetic CRS:  WGS 84
First 10 features:
   from_id  to_id travel_time_p10 travel_time_p50 parent_station_name                  geometry
1   100060 100060               0               0   Josef-Lang-Straße POINT (11.47202 48.15011)
2   100060  10157              28              NA   Josef-Lang-Straße POINT (11.46016 48.18474)
3   100060 105596              29              NA   Josef-Lang-Straße POINT (11.48913 48.18399)
4   100060 109574              23              NA   Josef-Lang-Straße POINT (11.47472 48.12805)
5   100060  11520              25              NA   Josef-Lang-Straße POINT (11.54544 48.14428)
6   100060 128167              23              NA   Josef-Lang-Straße POINT (11.46642 48.19177)
7   100060 131549              21              25   Josef-Lang-Straße POINT (11.44335 48.16361)
8   100060 131638              22              NA   Josef-Lang-Straße POINT (11.47278 48.17485)
9   100060 132365              19              25   Josef-Lang-Straße  POINT (11.4531 48.14563)
10  100060   1335              29              NA   Josef-Lang-Straße POINT (11.45716 48.20257)

This data format is easier to interpret if we look at a single stop. For example, lets find all stops that can be reached from Theresienstraße within a median travel time of <= 20 minutes.

ts_stops_20 <- filter(TTM,parent_station_name == "Theresienstraße", travel_time_p50 <= 20)


ggplot()+
  geom_sf(data = area)+
  geom_sf(data = ts_stops_20, color = "lightblue")+
  geom_sf(data = filter(stops,parent_station_name == "Theresienstraße"), color = "blue")

Let’s export this (& the stops) so that we can work with it in QGIS.

st_write(TTM,"./output/TTM.gpkg",append = F)

Deleting layer `TTM' using driver `GPKG'
Writing layer `TTM' to data source `./output/TTM.gpkg' using driver `GPKG'
Writing 295492 features with 5 fields and geometry type Point.

st_write(stops, "./output/stops.gpkg", append = F)

Deleting layer `stops' using driver `GPKG'
Writing layer `stops' to data source `./output/stops.gpkg' using driver `GPKG'
Writing 1125 features with 5 fields and geometry type Point.

We can achieve the same result as above by doing the following:

Then, by using what you’ve learned in previous tutorials, you can create service areas:

We’ve just scratched the surface of what you can do with r5r and tidytransit. Consider taking a look through the documentation for some inspiration for your project!